###############################################################################################
######### Script producing the tables and figures for the main analysis #######################
######### Bureaucratic Turnover Under New Governments: ########################################
######### The moderating effect of organizational buffering ###################################
###############################################################################################
####preamble####
#packages
library(dplyr)
library(survival)
library(survminer)
library(stargazer)
library(lmtest)
library(coxed)
library(ggthemes)
library(viridis)
library(ggpubr)
library(splitstackshape)
library(reshape2)
library(car)
library(sandwich)
library(lubridate)
library(stringr)
library(fixest)
library(marginaleffects)
library(tidyr)
library(modelsummary)

#clean
rm(list = ls())

#setting directory for output
wd_out <- "figures/main/"

#### functions ####
# creating figures
print_plot <- function(plot=obj, name, width=12, height=8){
  ggsave(paste0(name, ".pdf"), plot = plot, width = width, height = height, 
         units = "in", path = paste0("./", wd_out))
  ggsave(paste0(name, ".png"), plot = plot, width = width, height = height, 
         units = "in", path = paste0("./", wd_out,"pngs"))
}

# creating tex tables
output <- function(name, format = "tex") {
  paste0(wd_out, name, ".", format)
}

#Function to calculate marginal effects for the stepwise models with 1 interaction term between change of government and a dummy variable
model_marginal_effects <- function(list, NO_FE = 0) {
  for (i in 1:length(list)) {
    for(conf in c(0.9, 0.95)) {
      temp <- marginaleffects::avg_slopes(list[[i]], 
                                          variables = names(list[[i]]$coefficients)[1+NO_FE], 
                                          by = names(list[[i]]$coefficients)[2+NO_FE],
                                          conf_level = conf) %>% 
        mutate(!!names(.[3]) := case_when(names(.[3]) == "bureaucratic_buffer" & .[3] == 1 ~ "Level-2", 
                                          names(.[3]) == "bureaucratic_buffer" & .[3] == 0 ~ "Level-1",
                                          names(.[3]) == "state_secretaries_NSDID" & .[3] == 1 ~ "Level-1 & Political buffer",
                                          names(.[3]) == "state_secretaries_NSDID" & .[3] == 0 ~ "Level-1 & No buffer",
                                          TRUE ~ NA),
               c_lvl = conf)
      difference <- marginaleffects::avg_slopes(list[[i]], 
                                                variables = names(list[[i]]$coefficients)[1+NO_FE], 
                                                by = names(list[[i]]$coefficients)[2+NO_FE],
                                                hypothesis = "sequential",
                                                conf_level = conf) %>% 
        select(p.value, estimate, std.error, conf.low, conf.high) %>% 
        mutate(!!names(list[[1]]$coefficients)[2+NO_FE] := "Difference-in-Differences",
               term = names(list[[i]]$coefficients)[1+NO_FE],
               c_lvl = conf)
      
      temp <- bind_rows(temp, difference) %>% 
        mutate(model = names(list[i]))
      
      if (i == 1 & conf == 0.9) {
        all_models <- temp
      } else {
        all_models <- full_join(all_models, temp)
      }
    }
  }
  return(all_models)
  
}

#Function to calculate marginal effects for more complex models, specified hypotheses and conditional variable
model_marginal_effects_mechanisms <- function(model, conditional_variables_str, hyp) {
  for(conf in c(0.9, 0.95)) {
    temp <- marginaleffects::avg_slopes(model, 
                                        variables = names(model$coefficients)[1], 
                                        by = conditional_variables_str,
                                        conf_level = conf) %>% 
      mutate(!!names(.[3]) := case_when(names(.[3]) == "bureaucratic_buffer" & .[3] == 1 ~ "Level-2", 
                                        names(.[3]) == "bureaucratic_buffer" & .[3] == 0 ~ "Level-1",
                                        names(.[3]) == "state_secretaries_NSDID" & .[3] == 1 ~ "Level-1 & Political buffer",
                                        names(.[3]) == "state_secretaries_NSDID" & .[3] == 0 ~ "Level-1 & No buffer", 
                                        names(.[3]) == "election_year" & .[3] == 1 ~ "Election year",
                                        names(.[3]) == "election_year" & .[3] == 0 ~ "Within-term",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-1 (Permanent Secretary)" ~ "Level-1 (Permanent Secretary)",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-2 (Director General)" ~ "Level-2 (Director General)",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-1 (Director General)" ~ "Level-1 (Director General)",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-1 NO POL" ~ "Level-1 & No Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-2 and POL" ~ "Level-2 & Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-1 and POL" ~ "Level-1 & Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-2 NO POL" ~ "Level-2 & No Political buffer",
                                        .default = NA),
             c_lvl = conf)
    if (names(temp)[4] != "estimate") { 
      temp <- temp %>% 
        mutate(!!names(.[4]) := case_when(names(.[4]) == "bureaucratic_buffer" & .[4] == 1 ~ "Level-2", 
                                          names(.[4]) == "bureaucratic_buffer" & .[4] == 0 ~ "Level-1",
                                          names(.[4]) == "state_secretaries_NSDID" & .[4] == 1 ~ "Level-1 & Political buffer",
                                          names(.[4]) == "state_secretaries_NSDID" & .[4] == 0 ~ "Level-1 & No buffer", 
                                          names(.[4]) == "election_year" & .[4] == 1 ~ "Election year",
                                          names(.[4]) == "election_year" & .[4] == 0 ~ "Within-term",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-1 (Permanent Secretary)" ~ "Level-1 (Permanent Secretary)",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-2 (Director General)" ~ "Level-2 (Director General)",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-1 (Director General)" ~ "Level-1 (Director General)",
                                          .default = NA))
    }
    
    difference <- marginaleffects::avg_slopes(model, 
                                              variables = names(model$coefficients)[1], 
                                              by = conditional_variables_str,
                                              hypothesis = hyp,
                                              conf_level = conf) %>% 
      mutate(!!names(temp)[3] := "Difference-in-Differences",
             c_lvl = conf)
    
    if ("election_year" %in% conditional_variables_str) { 
      difference <- difference %>% 
        mutate(term = purrr::reduce2(c("1\\)", "0\\)"), c("Election year)", "Within-term)"), .init = term, str_replace_all))
    } else if ("state_secretaries_NSDID"  %in% conditional_variables_str) {
      difference <- difference %>% 
        mutate(term = purrr::reduce2(c("1\\)", "0\\)"), c("Political buffer)", "No Political buffer)"), .init = term, str_replace_all))
    }
    
    temp <- bind_rows(temp, difference)
    
    if (conf == 0.9) {
      all_models <- temp
    } else {
      all_models <- full_join(all_models, temp)
    }
  }
  return(all_models)
}

#Marginal effects function for Cox and Logistic models
cox_glm_marginal_effects_function <-  function(model, moderating_variabel, mod) { 
  
  combcoef <- c(1, 1, 0)*model$coefficients["GOV_turnover"]+
    c(0, 1, 1)*model$coefficients[str_c("GOV_turnover:", moderating_variabel)]
  
  if (mod == "Logistic") {
    vcov <- vcov(model, cluster ~ PersonID)
  } else { 
    vcov <- model$var 
    rownames(vcov) <- names(coefficients(model))
    colnames(vcov) <- names(coefficients(model))
  }
  
  combcoef.se <- c(sqrt(vcov["GOV_turnover","GOV_turnover"]),
                   sqrt(vcov["GOV_turnover","GOV_turnover"]+
                          (1^2)*vcov[str_c("GOV_turnover:", moderating_variabel), str_c("GOV_turnover:", moderating_variabel)]+
                          2*vcov["GOV_turnover", str_c("GOV_turnover:", moderating_variabel)]),
                   sqrt(vcov[str_c("GOV_turnover:", moderating_variabel), str_c("GOV_turnover:", moderating_variabel)]))
  
  group <- if (moderating_variabel == "bureaucratic_buffer") {
    group <- c("Level-1", "Level-2", "Difference-in-Differences")
  } else {
    group <- c("Level-1 & No buffer", "Level-1 & Political buffer", "Difference-in-Differences")
  }
  
  plotdata <- data.frame(name=moderating_variabel,
                         value=group,
                         estimate=exp(combcoef), 
                         conf.low=c(exp(combcoef+1.645*combcoef.se), exp(combcoef+1.96*combcoef.se)),
                         conf.high=c(exp(combcoef-1.645*combcoef.se), exp(combcoef-1.96*combcoef.se)),
                         c_lvl=c(rep("0.9", 3), rep("0.95", 3)),
                         term="Change of government",
                         model=mod)
  
  return(plotdata)
  
}

#Function to plot the marginal effects and the differences
marginal_effect_plot <- function(dataframe, facet = TRUE, models = FALSE, ORHR = FALSE, base_size = 16) {
  if (is.null(dataframe$name)) {
    dataframe <- dataframe %>% 
      pivot_longer(!!names(.[3]), values_drop_na = T)
  }
  
  if (facet == TRUE & models == FALSE & ORHR == FALSE) {
    plot <- ggplot(dataframe, aes(estimate, str_replace_all(term, output_names), xmin = conf.low, xmax = conf.high, 
                                  col = factor(value), linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = value))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(0.4, 1.2))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else if (facet == TRUE & models == TRUE & ORHR == FALSE) {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high, 
                                  col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = model))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(0.4, 1.2))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else if (facet == TRUE & models == TRUE & ORHR == TRUE) {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high, 
                                  col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = model))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 1, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect (Odds/Hazard Ratio)", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(1.2, 2.2))+
      scale_x_continuous(breaks = seq(0,6, by=0.5))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high,
                                  group = factor(value),
                                  linewidth = factor(c_lvl, levels = c("0.95", "0.9"))))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 3) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(1.4, 2.4))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  }
}

#### Lists of names ####
output_names <- c("GOV_turnover_lead" = "Change of government (t+1)",
                  "GOV_turnover_lagged" = "Change of government (t-1)",
                  "GOV_turnover" = "Change of government",
                  "minister_turnover" = "Change of minister",
                  "bureaucratic_buffer_i" = "Bureaucratic buffer & position type",
                  "bureaucratic_buffer" = "Bureaucratic buffer",
                  "organizational_buffer2" = "Organizational buffer",
                  "state_secretaries_NSDID" = "Political buffer",
                  "state_secretaries" = "Political buffer",
                  "years_since_last_turnover_2" = "Time since change of government",
                  "years_since_last_turnover_2_categories" = "Time since change of government categorial",
                  "match_start_end_PMparty_flipped" = "Party incongruence",
                  "election_year" = "Election year",
                  "age_year" = "Age",
                  "gender" = "Gender",
                  "temporary_employment" = "Temporary appointment",
                  "pre_unionsoppløsning" = "Before-1906", 
                  "foreign_affairs" = "Ministry of Foreign Affairs",
                  "government_PMparty_share_posts_end_of_year" = "PM-party share of cabinet",
                  "years_since_DEP_start" = "Time since ministry created",
                  "dep_end" = "Ministry terminated",
                  "pol_flexible" = "Exposed to change of government",
                  "NSD_ID" = "Ministry Fixed Effects",
                  "decade" = "Decade",
                  "decade2" = "Decade",
                  "duration" = "Senior bureaucrat tenure in years",
                  "WWII" = "WWII",
                  "age_group_year2" = "Age Group",
                  "multiple_L1" = "Multiple Level-1 bureuacrats",
                  "dep.rad_turnover" = "Senior bureaucrat turnover")

moderators_names <- str_replace_all(sapply(c("Baseline",
                                             "government_PMparty_share_posts_end_of_year",
                                             "years_since_DEP_start", "dep_end",
                                             "multiple_L1", "age_group_year2"), 
                                           function(x){rep(x, 6)}), output_names)

DV_names <- lapply(c("Senior bureaucrat turnover", 
                     "Including ministry switching",
                     "Including ministry switching \n and promotion/demotion"), 
                   function(x){rep(x, 6)}) %>% unlist()

#### modelsummary fit statistics ####
gm <- tribble(~raw, ~clean, ~fmt,
              "FE: time", "FE: Time", 0,
              "FE: NSD_ID", "FE: Ministry", 0,
              "FE: decade", "FE: Decade", 0,
              "nobs", "N", 0,
              "r.squared", "R2", 2,
              "aic", "AIC", 2)

###################################
####preparing data for analysis####
###################################
#Load data
df <- readRDS(file = "data/TCS_NOR.RDS") 

#recoding
df$match_start_end_PMparty_flipped <- if_else(df$match_start_end_PMparty == 1, 0, 1)
df$election_year <- if_else(is.na(df$election_year), 0, df$election_year)
df$military <- ifelse(is.na(df$military), 0, df$military)



#### creating subsets of the data ####
raw <- df %>% filter(year >= 1884, military == 0) #removing military DGs and observations before parliamentarism

#pensioners are removed | and WWII
df <- raw %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70, WWII == 0) 


#pensioners are removed 
df_med_WWII <- raw %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70) 

#spells must be more than 2 years and pensioners are removed 
df2 <- df %>% 
  filter(less_than_two_years == 0) %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70) 

#only Level-1 bureaucrats
df_level1 <- df %>% filter(bureaucratic_buffer == 0)

#### dataframe with all used variables
sub_df <- df_med_WWII %>% 
  select(start, stop, dep.rad_turnover, GOV_turnover, years_since_last_turnover_2, match_start_end_PMparty_flipped,
         election_year, age_year, gender, temporary_employment, pol_flexible, pre_unionsoppløsning, WWII, decade2,
         dep_end, years_since_DEP_start, government_PMparty_share_posts_end_of_year,
         foreign_affairs, WWI, post_WWII, NSD_ID, bureaucratic_buffer, year, 
         state_secretaries, proximity_to_GOV_turnover,proximity_to_GOV_turnover_2_categories, years_since_last_turnover_2_categories,
         GOV_turnover_any_change, GOV_turnover_new_PM, GOV_turnover_new_Parties,
         GOV_turnover_during_term, GOV_turnover_election, government_percent_parliament_end_of_year, 
         government_single_party_end_of_year, government_minority_end_of_year, government_left_right_gov_end_of_year,
         gov_turnover_during_spell, pension_age_propensity_age_year, pension_age_propensity_age_at_start,
         match_start_end_party, match_start_end_left_right, match_start_end_coalition_party,
         minister_turnover, political_career,
         GOV_turnover_lagged, GOV_turnover_lead, id_spells, duration, bureaucratic_buffer_i, 
         y_since_PS_in_min, decade_since_PS_in_min,rounded_y_since_PS_in_min, organizational_buffer, 
         organizational_buffer2, after_political_appointees_NSDID, after_state_secretaries_NSDID, after_poladv_NSDID,
         next_position_sector, education_main, education_level, post_SLS, pension_age_propensity_age_year,
         age_group_year2, multiple_L1, ministry_portfolio)

#### Statistics mentioned in text ####
sum(df$age_year == 70) # number of senior bureaucrats who reached mandatory retirment
length(unique(df_med_WWII$PersonID)) #number of people
length(unique(df_med_WWII$id_spells)) #number of employment spells 
sum(df_med_WWII$dep.rad_turnover == 1) #number of employment spells that ended in the study period
length(unique(df_med_WWII$PersonID[df_med_WWII$bureaucratic_buffer == 0])) #number of people serving as level-1 bureaucrats
length(unique(df_med_WWII$PersonID[df_med_WWII$n_spells > 1])) #number of people with more than one employment spells

#### Figure 1 - The number of senior bureaucrats in Norwegian ministries 1884-2024 subset by position####
figure1 <- ggplot(df_med_WWII %>%
                      distinct(id_spells, year, bureaucratic_buffer_i), 
                    aes(year, fill = factor(bureaucratic_buffer_i)))+
  geom_bar(width = 0.9)+
  scale_x_continuous(limits = c(1883, 2025), breaks = seq(1884, 2024, 5))+
  scale_y_continuous(limits = c(0, 150))+
  theme_base(base_size = 14)+
  scale_fill_viridis(discrete = T,
                     breaks = c("Level-1 (Permanent Secretary)",
                                "Level-1 (Director General)",
                                "Level-2 (Director General)"),
                     direction = -1)+
  labs(fill = NULL, y = "Number of Senior Bureaucrats", x = NULL) +
  theme(legend.position = "top",
        legend.box = "vertical", 
        plot.background=element_blank(), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

print_plot(plot = figure1, name = "figure1")


#### Figure 2 - The proportion of senior bureaucratic turnover 1884-2024 ####
figure2 <- ggplot(df_med_WWII %>%
                      distinct(year, mean.dep.rad.turnoverrate.year, GOV_turnover) %>% 
                    mutate(mean_total = mean(mean.dep.rad.turnoverrate.year)),
                    aes(x = year, y = mean.dep.rad.turnoverrate.year, fill = as.factor(GOV_turnover)))+
  geom_col()+
  scale_x_continuous(limits = c(1883, 2025), breaks = seq(1884, 2024, 5))+
  scale_y_continuous(limits = c(0, 1), breaks = c(0, .11, .25, .5, .75, 1))+
  scale_fill_viridis(discrete = T, labels = c("No", "Yes"), breaks = c(0, 1))+
  geom_hline(aes(yintercept=mean_total), linetype=2)+
  theme_base(base_size = 16)+
  labs(fill = "Change of Government", y = "Proportion of Exits", x = NULL) +
  theme(legend.position = "top",
        legend.box = "vertical", 
        plot.background=element_blank(), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

print_plot(plot = figure2, name = "figure2")


#### Figure 3 - Boxplot of the proportion of senior bureaucrat turnover in years with and without change of government subset by organizational buffering ####
bp_data <- full_join(df %>% 
                       filter(WWII == 0) %>% 
                       distinct(year, id_spells, dep.rad_turnover, GOV_turnover, bureaucratic_buffer) %>% 
                       summarise(mean_turnover = mean(dep.rad_turnover), 
                                 .by = c(year, bureaucratic_buffer, GOV_turnover)) %>% 
                       mutate(x_variabel = if_else(bureaucratic_buffer == 1, "Level-2", "Level-1"),
                              grp_variabel = "Bureaucratic buffer") %>% 
                       select(-bureaucratic_buffer),
                     df %>% 
                       filter(WWII == 0, bureaucratic_buffer == 0) %>% 
                       distinct(year, id_spells, dep.rad_turnover, GOV_turnover, state_secretaries) %>% 
                       summarise(mean_turnover = mean(dep.rad_turnover), 
                                 .by = c(year, state_secretaries, GOV_turnover)) %>% 
                       mutate(x_variabel = if_else(state_secretaries == 1, "Level-1 & Political buffer", "Level-1 & No buffer"),
                              grp_variabel = "Political buffer") %>% 
                       select(-state_secretaries))

#to write about the descriptives
bp_data %>% 
  distinct(year,mean_turnover,
           GOV_turnover, x_variabel, grp_variabel) %>% 
  summarise(mean_turnover = mean(mean_turnover), .by = c(year, x_variabel, grp_variabel, GOV_turnover)) %>% 
  summarise(mean = mean(mean_turnover), .by = c(x_variabel, grp_variabel, GOV_turnover))

figure3 <- ggplot(bp_data,
                   aes(x = as.factor(x_variabel), y = mean_turnover, fill = as.factor(GOV_turnover)))+
  geom_boxplot()+
  scale_x_discrete()+
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.10))+
  scale_fill_viridis(discrete = T, labels = c("No", "Yes"), direction = 1)+
  theme_base(base_size = 16)+
  labs(fill = "Change of Government", y = "Proportion of Exits", x = NULL) +
  theme(legend.position = "top",
        legend.box = "vertical",
        plot.background=element_blank())+
  facet_wrap(~grp_variabel, drop = T, scales = "free_x")

print_plot(plot = figure3, name = "figure3")

#### Table 1 ####
# Bureaucratic buffer main model
OLS_FE <- feols(dep.rad_turnover ~ GOV_turnover * bureaucratic_buffer +
                  csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | time+NSD_ID+decade, 
                cluster = "PersonID", df)
modelsummary(OLS_FE, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_FE", "html"), 
             coef_rename = output_names,
             col.names = NULL,
             gof_map = gm,
             title = "Step-wise OLS regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

# Political buffer main model
OLS_FE_polbuff <- feols(dep.rad_turnover ~ GOV_turnover * state_secretaries_NSDID +
                          csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | time+NSD_ID+decade, 
                        cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
modelsummary(OLS_FE_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_FE_polbuff", "html"), 
             coef_rename = output_names,
             col.names = NULL,
             gof_map = gm,
             title = "Step-wise OLS regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Calculating marginal effects and diff-in-diff
Coefplot_OLS_FE <- model_marginal_effects(OLS_FE)
Coefplot_OLS_FE_polbuff <- model_marginal_effects(OLS_FE_polbuff)

main_table <- bind_rows(
  tibble(bureaucratic_buffer = "Panel 1: Bureaucratic Buffer", model = " ", 
         dvest = mean(df$dep.rad_turnover, na.rm = T),
         N = str_c(nrow(df), " [", length(unique(df$id_spells)), "] (", sum(df$dep.rad_turnover == 1), ")")
         ),
  Coefplot_OLS_FE %>% 
    filter(nchar(model) == max(nchar(model))) %>% 
    mutate(model = "OLS FE"),
  tibble(state_secretaries_NSDID = " ", model = " "),
  tibble(state_secretaries_NSDID = "Panel 2: Political Buffer", model = " ",
         dvest = mean(df$dep.rad_turnover[df$bureaucratic_buffer == 0], na.rm = T),
         N = str_c(nrow(df_level1), " [", length(unique(df_level1$id_spells)), "] (", sum(df_level1$dep.rad_turnover == 1), ")")
  ),
  Coefplot_OLS_FE_polbuff %>% 
    filter(nchar(model) == max(nchar(model))) %>% 
    mutate(model = "OLS FE") ) %>% 
  pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T, values_to = "variable") %>% 
  select(`Change of government (Marginal effect)` = variable, estimate, std.error, "Average senior bureaucrat turnover probability" = dvest, p.value, N) %>% 
  distinct() %>% 
  mutate(across(estimate:p.value, ~round(.x, digits = 3))) %>% 
  mutate(across(estimate:N, ~replace_na(as.character(.x), " ")))

notes_main <- c("Note: The table reports Marginal effects of change of government obtained from OLS models",
                "with fixed effects for tenure, decade and ministry. The model also includes control",
                "variables for election year and government independence, as well as a set of individual char-",
                "acteristics: Age, gender, and temporary appointment. The WWII years are dropped from all models.",
                "Standard errors are clustered on individuals. The dependent variable is turnover of individual",
                "bureaucrats. The N-column displays the number of individual year observations, with the number of",
                "employment spells in square brackets and the number of bureaucratic turnover events in parentheses.", 
                "Variation in N is due to panel 2 only including Level-1 bureaucrats.")

stargazer(main_table,
          summary = F,
          title = "Marginal effect of change of government from main OLS FE regression models",
          align=F,
          no.space=T, 
          model.numbers=T,
          rownames = F,
          label = c("Tab:main_marginal_table"),
          type = "latex",
          notes = notes_main,
          out = output(name = "table1_raw", format = "tex")
) #The layout of the table is further customized in LaTeX

#HMTL version of the two regression models used to create the main table
modelsummary(list(OLS_FE[[6]], OLS_FE_polbuff[[6]]), 
             conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("main_results", "html"), 
             coef_rename = output_names,
             coef_omit = c(-1,-2, -8, -9, -10),
             gof_map = gm,
             add_rows = data.frame(rbind(cbind("Number of senior bureaucrats", 
                                               length(unique(df$id_spells[is.finite(df$bureaucratic_buffer)])),
                                               length(unique(df$id_spells[is.finite(df$bureaucratic_buffer) & df$bureaucratic_buffer_i != "Level-2 (Director General)"]))),
                                         cbind("Number of events", sum(df$dep.rad_turnover == 1), sum(df_level1$dep.rad_turnover == 1)))),
             title = "Regression results:",
             notes = list("Note: DV = Turnover of senior bureaucrat. Estimates reported with 95 percent confidence intervals (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)



#### Figure 4 - Marginal effects of change of government on bureaucratic turnover conditional on bureaucratic buffering and position type ####
table(df$bureaucratic_buffer_i) #N observations
OLS_different_positions <- feols(dep.rad_turnover ~ GOV_turnover*bureaucratic_buffer_i + election_year +
                                   age_year + gender + temporary_employment + pre_unionsoppløsning | NSD_ID + time + decade, 
                                 data = df)
figure4 <- model_marginal_effects_mechanisms(OLS_different_positions,
                                         conditional_variables_str = c("bureaucratic_buffer_i"),
                                         hyp = "revpairwise")
print_plot(marginal_effect_plot(figure4)+
             theme(strip.background = element_blank(),
                   strip.text.x = element_blank())+
             guides(color = guide_legend(reverse=TRUE, nrow = 2, position = "top")), 
           name = "figure4")

#### Figure 5 - Event study estimates of the changing marginal effect of change of government on bureaucratic turnover in 3 periods before and after the introduction of permanent secretaries as an organizational buffer in a ministry ####
table(df$rounded_y_since_PS_in_min, df$GOV_turnover, df$bureaucratic_buffer) #descriptive at risk at the different periods
table(df$rounded_y_since_PS_in_min, df$GOV_turnover, df$bureaucratic_buffer_i) #descriptive at risk at the different periods
table(df$bureaucratic_buffer_i) #descriptive at risk at the different periods
table(df$rounded_y_since_PS_in_min, df$y_since_PS_in_min)
OLS_event_gov_turnover <- feols(dep.rad_turnover ~ GOV_turnover * i(rounded_y_since_PS_in_min, ref=0, bureaucratic_buffer_i) + election_year +
                                  age_year + gender + temporary_employment + pre_unionsoppløsning | NSD_ID + time + decade, 
                                data = df %>% filter())
summary(OLS_event_gov_turnover)

pdata <- bind_rows(avg_slopes(OLS_event_gov_turnover, 
                              variables = "GOV_turnover", 
                              by = c("rounded_y_since_PS_in_min", "bureaucratic_buffer_i"), conf_level = 0.95) %>% 
                     mutate(c_lvl = "0.95", model = "B"),
                   avg_slopes(OLS_event_gov_turnover, 
                              variables = "GOV_turnover", 
                              by = c("rounded_y_since_PS_in_min", "bureaucratic_buffer_i"), conf_level = 0.90) %>% 
                     mutate(c_lvl = "0.90", model = "A"))

figure5 <- ggplot(pdata, aes(as.factor(rounded_y_since_PS_in_min), estimate, ymin = conf.low, ymax = conf.high, 
                         col = bureaucratic_buffer_i, linewidth = factor(c_lvl, levels = c("0.95", "0.90")),
                         group = bureaucratic_buffer_i))+
  geom_linerange(position = position_dodge(0.4))+
  geom_point(position = position_dodge(0.4), size = 3)+
  geom_hline(yintercept = 0, linetype = 2, color = "darkred")+
  geom_vline(xintercept = 3.5, linetype = 1, color = "black")+
  scale_x_discrete(labels = c("[min, -40)", "[-40,-20)", "[-20, 0)", "[0, 15]", "(15, 30]", "(30, max]"))+
  scale_linewidth_discrete(range = c(1.4, 2.4))+
  scale_y_continuous(breaks = seq(-0.2,0.2, by=0.05))+
  theme_base(base_size = 18)+
  scale_color_viridis(discrete = T, name = "", option = "D")+
  labs(x = "Years since introduction of permanent secretary as level-1", y = "Marginal effect of change of government") +
  theme(legend.position = "top",
        legend.box = "vertical", plot.background=element_blank())+
  guides(color = guide_legend(reverse=TRUE),
         linewidth = "none")

print_plot(plot = figure5, name = "figure5")


#### Figure 6 - Marginal effects of change of government on bureaucratic turnover conditional on bureaucratic buffering and elections ####
OLS_election_year <- feols(dep.rad_turnover ~ GOV_turnover*election_year*bureaucratic_buffer +
                             age_year + gender + temporary_employment + pre_unionsoppløsning | NSD_ID + time + decade, 
                           data = df)
obj <- model_marginal_effects_mechanisms(OLS_election_year,
                                         conditional_variables_str = c("bureaucratic_buffer", "election_year"),
                                         hyp = "revpairwise")

# ploting marginal effects compared to Within-term
figure6 <- obj %>% 
  mutate(`bureaucratic_buffer & state_secretaries_NSDID` = str_remove_all(paste(bureaucratic_buffer, election_year, sep = " & "), " & NA"),
         .after = "contrast") %>% 
  filter(term == "GOV_turnover")
print_plot(marginal_effect_plot(figure6, facet = F), name = "figure6")

#### Figure 7 - Marginal effects with 95\% (small stroke width) and 90\% (big stroke width) confidence bands for the marginal effect of different operationalizations of change of government ####
#bureaucratic buffering
OLS_alternative_IVs <- feols(dep.rad_turnover ~ sw(GOV_turnover * bureaucratic_buffer, 
                                                   GOV_turnover_lagged * bureaucratic_buffer, 
                                                   GOV_turnover_lead * bureaucratic_buffer,
                                                   minister_turnover * bureaucratic_buffer,
                                                   match_start_end_PMparty_flipped * bureaucratic_buffer) +
                               election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade, 
                             cluster = "PersonID", df)
names(OLS_alternative_IVs) <- names(OLS_alternative_IVs) <- str_replace_all(c("GOV_turnover", "GOV_turnover_lagged", "GOV_turnover_lead", "minister_turnover", "match_start_end_PMparty_flipped"), output_names)
modelsummary(OLS_alternative_IVs, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_alternative_IVs", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: Different IV operationalizations",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

#political buffering
OLS_alternative_IVs_polbuff <- feols(dep.rad_turnover ~ sw(GOV_turnover * state_secretaries_NSDID, 
                                                           GOV_turnover_lagged * state_secretaries_NSDID, 
                                                           GOV_turnover_lead * state_secretaries_NSDID,
                                                           minister_turnover * state_secretaries_NSDID,
                                                           match_start_end_PMparty_flipped * state_secretaries_NSDID) +
                                       election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade, 
                                     cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
names(OLS_alternative_IVs_polbuff) <- names(OLS_alternative_IVs_polbuff) <- str_replace_all(c("GOV_turnover", "GOV_turnover_lagged", "GOV_turnover_lead", "minister_turnover", "match_start_end_PMparty_flipped"), output_names)
modelsummary(OLS_alternative_IVs_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_alternative_IVs_polbuff", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: Different IV operationalizations",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

#calculating marginal effects and outputing
figure7 <- model_marginal_effects(OLS_alternative_IVs_polbuff) %>% 
  full_join(., model_marginal_effects(OLS_alternative_IVs)) %>% 
  pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T)
print_plot(marginal_effect_plot(figure7), name = "figure7")
